home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zhetrd.f < prev    next >
Text File  |  1997-06-25  |  9KB  |  282 lines
  1.       SUBROUTINE ZHETRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          UPLO
  10.       INTEGER            INFO, LDA, LWORK, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   D( * ), E( * )
  14.       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  ZHETRD reduces a complex Hermitian matrix A to real symmetric
  21. *  tridiagonal form T by a unitary similarity transformation:
  22. *  Q**H * A * Q = T.
  23. *
  24. *  Arguments
  25. *  =========
  26. *
  27. *  UPLO    (input) CHARACTER*1
  28. *          = 'U':  Upper triangle of A is stored;
  29. *          = 'L':  Lower triangle of A is stored.
  30. *
  31. *  N       (input) INTEGER
  32. *          The order of the matrix A.  N >= 0.
  33. *
  34. *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  35. *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
  36. *          N-by-N upper triangular part of A contains the upper
  37. *          triangular part of the matrix A, and the strictly lower
  38. *          triangular part of A is not referenced.  If UPLO = 'L', the
  39. *          leading N-by-N lower triangular part of A contains the lower
  40. *          triangular part of the matrix A, and the strictly upper
  41. *          triangular part of A is not referenced.
  42. *          On exit, if UPLO = 'U', the diagonal and first superdiagonal
  43. *          of A are overwritten by the corresponding elements of the
  44. *          tridiagonal matrix T, and the elements above the first
  45. *          superdiagonal, with the array TAU, represent the unitary
  46. *          matrix Q as a product of elementary reflectors; if UPLO
  47. *          = 'L', the diagonal and first subdiagonal of A are over-
  48. *          written by the corresponding elements of the tridiagonal
  49. *          matrix T, and the elements below the first subdiagonal, with
  50. *          the array TAU, represent the unitary matrix Q as a product
  51. *          of elementary reflectors. See Further Details.
  52. *
  53. *  LDA     (input) INTEGER
  54. *          The leading dimension of the array A.  LDA >= max(1,N).
  55. *
  56. *  D       (output) DOUBLE PRECISION array, dimension (N)
  57. *          The diagonal elements of the tridiagonal matrix T:
  58. *          D(i) = A(i,i).
  59. *
  60. *  E       (output) DOUBLE PRECISION array, dimension (N-1)
  61. *          The off-diagonal elements of the tridiagonal matrix T:
  62. *          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
  63. *
  64. *  TAU     (output) COMPLEX*16 array, dimension (N-1)
  65. *          The scalar factors of the elementary reflectors (see Further
  66. *          Details).
  67. *
  68. *  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
  69. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  70. *
  71. *  LWORK   (input) INTEGER
  72. *          The dimension of the array WORK.  LWORK >= 1.
  73. *          For optimum performance LWORK >= N*NB, where NB is the
  74. *          optimal blocksize.
  75. *
  76. *  INFO    (output) INTEGER
  77. *          = 0:  successful exit
  78. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  79. *
  80. *  Further Details
  81. *  ===============
  82. *
  83. *  If UPLO = 'U', the matrix Q is represented as a product of elementary
  84. *  reflectors
  85. *
  86. *     Q = H(n-1) . . . H(2) H(1).
  87. *
  88. *  Each H(i) has the form
  89. *
  90. *     H(i) = I - tau * v * v'
  91. *
  92. *  where tau is a complex scalar, and v is a complex vector with
  93. *  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
  94. *  A(1:i-1,i+1), and tau in TAU(i).
  95. *
  96. *  If UPLO = 'L', the matrix Q is represented as a product of elementary
  97. *  reflectors
  98. *
  99. *     Q = H(1) H(2) . . . H(n-1).
  100. *
  101. *  Each H(i) has the form
  102. *
  103. *     H(i) = I - tau * v * v'
  104. *
  105. *  where tau is a complex scalar, and v is a complex vector with
  106. *  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
  107. *  and tau in TAU(i).
  108. *
  109. *  The contents of A on exit are illustrated by the following examples
  110. *  with n = 5:
  111. *
  112. *  if UPLO = 'U':                       if UPLO = 'L':
  113. *
  114. *    (  d   e   v2  v3  v4 )              (  d                  )
  115. *    (      d   e   v3  v4 )              (  e   d              )
  116. *    (          d   e   v4 )              (  v1  e   d          )
  117. *    (              d   e  )              (  v1  v2  e   d      )
  118. *    (                  d  )              (  v1  v2  v3  e   d  )
  119. *
  120. *  where d and e denote diagonal and off-diagonal elements of T, and vi
  121. *  denotes an element of the vector defining H(i).
  122. *
  123. *  =====================================================================
  124. *
  125. *     .. Parameters ..
  126.       DOUBLE PRECISION   ONE
  127.       PARAMETER          ( ONE = 1.0D+0 )
  128.       COMPLEX*16         CONE
  129.       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
  130. *     ..
  131. *     .. Local Scalars ..
  132.       LOGICAL            UPPER
  133.       INTEGER            I, IINFO, IWS, J, KK, LDWORK, NB, NBMIN, NX
  134. *     ..
  135. *     .. External Subroutines ..
  136.       EXTERNAL           XERBLA, ZHER2K, ZHETD2, ZLATRD
  137. *     ..
  138. *     .. Intrinsic Functions ..
  139.       INTRINSIC          MAX
  140. *     ..
  141. *     .. External Functions ..
  142.       LOGICAL            LSAME
  143.       INTEGER            ILAENV
  144.       EXTERNAL           LSAME, ILAENV
  145. *     ..
  146. *     .. Executable Statements ..
  147. *
  148. *     Test the input parameters
  149. *
  150.       INFO = 0
  151.       UPPER = LSAME( UPLO, 'U' )
  152.       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  153.          INFO = -1
  154.       ELSE IF( N.LT.0 ) THEN
  155.          INFO = -2
  156.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  157.          INFO = -4
  158.       ELSE IF( LWORK.LT.1 ) THEN
  159.          INFO = -9
  160.       END IF
  161.       IF( INFO.NE.0 ) THEN
  162.          CALL XERBLA( 'ZHETRD', -INFO )
  163.          RETURN
  164.       END IF
  165. *
  166. *     Quick return if possible
  167. *
  168.       IF( N.EQ.0 ) THEN
  169.          WORK( 1 ) = 1
  170.          RETURN
  171.       END IF
  172. *
  173. *     Determine the block size.
  174. *
  175.       NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
  176.       NX = N
  177.       IWS = 1
  178.       IF( NB.GT.1 .AND. NB.LT.N ) THEN
  179. *
  180. *        Determine when to cross over from blocked to unblocked code
  181. *        (last block is always handled by unblocked code).
  182. *
  183.          NX = MAX( NB, ILAENV( 3, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
  184.          IF( NX.LT.N ) THEN
  185. *
  186. *           Determine if workspace is large enough for blocked code.
  187. *
  188.             LDWORK = N
  189.             IWS = LDWORK*NB
  190.             IF( LWORK.LT.IWS ) THEN
  191. *
  192. *              Not enough workspace to use optimal NB:  determine the
  193. *              minimum value of NB, and reduce NB or force use of
  194. *              unblocked code by setting NX = N.
  195. *
  196.                NB = MAX( LWORK / LDWORK, 1 )
  197.                NBMIN = ILAENV( 2, 'ZHETRD', UPLO, N, -1, -1, -1 )
  198.                IF( NB.LT.NBMIN )
  199.      $            NX = N
  200.             END IF
  201.          ELSE
  202.             NX = N
  203.          END IF
  204.       ELSE
  205.          NB = 1
  206.       END IF
  207. *
  208.       IF( UPPER ) THEN
  209. *
  210. *        Reduce the upper triangle of A.
  211. *        Columns 1:kk are handled by the unblocked method.
  212. *
  213.          KK = N - ( ( N-NX+NB-1 ) / NB )*NB
  214.          DO 20 I = N - NB + 1, KK + 1, -NB
  215. *
  216. *           Reduce columns i:i+nb-1 to tridiagonal form and form the
  217. *           matrix W which is needed to update the unreduced part of
  218. *           the matrix
  219. *
  220.             CALL ZLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK,
  221.      $                   LDWORK )
  222. *
  223. *           Update the unreduced submatrix A(1:i-1,1:i-1), using an
  224. *           update of the form:  A := A - V*W' - W*V'
  225. *
  226.             CALL ZHER2K( UPLO, 'No transpose', I-1, NB, -CONE,
  227.      $                   A( 1, I ), LDA, WORK, LDWORK, ONE, A, LDA )
  228. *
  229. *           Copy superdiagonal elements back into A, and diagonal
  230. *           elements into D
  231. *
  232.             DO 10 J = I, I + NB - 1
  233.                A( J-1, J ) = E( J-1 )
  234.                D( J ) = A( J, J )
  235.    10       CONTINUE
  236.    20    CONTINUE
  237. *
  238. *        Use unblocked code to reduce the last or only block
  239. *
  240.          CALL ZHETD2( UPLO, KK, A, LDA, D, E, TAU, IINFO )
  241.       ELSE
  242. *
  243. *        Reduce the lower triangle of A
  244. *
  245.          DO 40 I = 1, N - NX, NB
  246. *
  247. *           Reduce columns i:i+nb-1 to tridiagonal form and form the
  248. *           matrix W which is needed to update the unreduced part of
  249. *           the matrix
  250. *
  251.             CALL ZLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ),
  252.      $                   TAU( I ), WORK, LDWORK )
  253. *
  254. *           Update the unreduced submatrix A(i+nb:n,i+nb:n), using
  255. *           an update of the form:  A := A - V*W' - W*V'
  256. *
  257.             CALL ZHER2K( UPLO, 'No transpose', N-I-NB+1, NB, -CONE,
  258.      $                   A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE,
  259.      $                   A( I+NB, I+NB ), LDA )
  260. *
  261. *           Copy subdiagonal elements back into A, and diagonal
  262. *           elements into D
  263. *
  264.             DO 30 J = I, I + NB - 1
  265.                A( J+1, J ) = E( J )
  266.                D( J ) = A( J, J )
  267.    30       CONTINUE
  268.    40    CONTINUE
  269. *
  270. *        Use unblocked code to reduce the last or only block
  271. *
  272.          CALL ZHETD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ),
  273.      $                TAU( I ), IINFO )
  274.       END IF
  275. *
  276.       WORK( 1 ) = IWS
  277.       RETURN
  278. *
  279. *     End of ZHETRD
  280. *
  281.       END
  282.